function [BW,rawpow,filtpow] = cvc_erpower(EEG,t1,t2,NW,f0s,chs)
% CVC_ERPOWER
%
%  Synopsis
%  ========
%
%                  BW  = cvc_erpower(EEG,t1,t2,NW)
%  [BW,rawpow,filtpow] = cvc_erpower(EEG,t1,t2,NW,f0s)
%  [BW,rawpow,filtpow] = cvc_erpower(EEG,t1,t2,NW,f0s,chs)
%
%  -- Author: Mads Dyrholm --
%     Center for Visual Cognition, University of Copenhagen.
%     2010 - March 2011
%
%  Purpose
%  =======
%
%  Time-window and frequency-band power estimation.
%
%  Inputs
%  ======
%
%  EEG - EEGLAB dataset to analyze. Each epoch will be
%  analyzed separately. The number of epochs is N.
%
%  t1 - Analysis window start time in milliseconds.
%
%  t2 - Analysis window end time in milliseconds.
%
%  NW - For multi-taper spectral analysis, typical choices
%  for NW are 2, 5/2, 3, 7/2, or 4, see also DPSS.
%
%  f0s - Vector of length L with the center frequencies
%  to analyze. A zero indicates no modulation (lowpass).
%
%  chs - Vector of length M with the list of channels to analyze.
%  Set element negative (-#) for analysis of ICA component (#).
%  This argument is optional, if omitted then all channels will 
%  be analyzed.
%  
%  Outputs
%  =======
%
%  BW - bandwidth.
%
%  rawpow - M-by-N matrix with the raw data power in each channel
%  and epoch. (e.g. channels-by-epochs-by-frequencies)
%
%  filtpow - M-by-N-by-L matrix with the power in each channel, 
%  epoch, and frequency band.
%
%  fmids - Analysis center frequencies.
%
%  References
%  ==========
%
%  Dyrholm, M. et al. (2011). Single-Trial Inference 
%  on Visual Attention, in proceedings of Computational Models
%  for Life Sciences (CMLS-11), Tokyo, Japan
%
%  Slepian, D. (1978). The Bell Systems Tech. J. 57, 
%  1371-1430.
%
%  Thomson, D. J. (1982). Proc. IEEE 70, 1055-1096.

idx_base = find(EEG.times>t1 & EEG.times<t2);

%f1 = f1/EEG.srate;
%f2 = f2/EEG.srate;
N = length(idx_base);
[E,V] = dpss(N,NW);
M = size(E,2);
W = NW/N;
BW = W*EEG.srate;
fprintf('Window: %i samples, Slepian sequences: %i, bandwidth: %f Hz\n',N,M,BW);
if nargin<5
  return
end
if nargin<6
  chs = 1:EEG.nbchan;
end
fprintf('Analyzing bands:');
Win = [-BW/2,+BW/2];
fmids = (f0s/EEG.srate);
rawpow = [];
filtpow = [];
for fmid=fmids
  fprintf(' [%.2f,%.2f]',fmid*EEG.srate+Win(1),fmid*EEG.srate+Win(2));
  if fmid>0
    shift_vec = exp(i*2*pi*fmid*(0:N-1))';
    F = E.*repmat(shift_vec,1,M);
  else
    F = E;
  end
  
  filtpowf = [];
  for ch = chs
    if ch<0
      x = squeeze(EEG.icaact(-ch,:,:));
    else
      x = squeeze(EEG.data(ch,:,:));
    end
    
    x = double(x);
  
    if isempty(filtpow), rawpow = cat(1,rawpow,sum(x(idx_base,:).^2)); end
    
    pp = F' * x(idx_base,:);
    filtpowf = cat(1,filtpowf,real(sum(pp.*conj(pp))));
  end
  filtpow = cat(3,filtpow,filtpowf);
end
fprintf('\n');
